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ABSTRACT 

We calculate the meteorology of the close-in transiting extrasolar planet HD 209458b using a global, three- 
dimensional atmospheric circulation model. Dynamics are driven by perpetual irradiation of one hemisphere of 
this tidally locked planet. The simulation predicts global temperature contrasts of ~ 500 K at the photosphere 
and the development of a steady superrotating jet. The jet extends from the equator to mid-latitudes and from 
the top model layer at 1 mbar down to 10 bars at the base of the heated region. Wind velocities near the equator 
exceed 4 km s" 1 at 300 mbar. The hottest regions of the atmosphere are blown downstream from the substellar 
point by ~ 60° of longitude. We predict from these results a factor of ~ 2 ratio between the maximum and 
minimum observed radiation from the planet over a full orbital period, with peak infrared emission preceding 
the time of the secondary eclipse by ~ 14 hours. 

Subject headings: planets and satellites: general — planets and satellites: individual (HD 209458b) — methods: 
numerical — atmospheric effects 



1. INTRODUCTION 

The transiting planet HD 209458b orbits very closely 
(0.046 AU) to its parent star with a period of 3.5257 days 
jCharbonneau et al J 1200(1 iHenrv et al J 120001) . From transit 
depth measurements, the mass and radius of HD 209458b 
are known f airly accurately: 0. 69 ± 0.05 Mj U pj te r and 1.32 ± 
0.05/?j up itei (Laughlin et al. 2005). The age of the system is 
estimated to be 5.2 Gyr, with uncertainties of ~ 10%. Fur- 
thermore, owing to careful measurements of the stellar spec- 
trum during the planet's transit, much i s now known about 
■ HP 209 45 8b 's atmospheric properties Ipr own et all 12001: 
ICharbonneau et alJ2002t IVidal-Ma dlar et all2003l 120041) . 

Considerable work has been done to model the spectra, 
physical structur e, and time evolution of extrasolar giant 
planets (EGPs ) (Burrows eTaTl 120041 iChabrier et all [2004; 
llro et alJ2005l) . Relatively less effort, however, has been spent 
on EGP meteorologies; i.e., global temperature and pressure 
fluctuations, wind velocities, and cloud propertie s. Prelim- 
inary simulations of the circulation by Sho wman & Gu illot 
(12002) and the shallow-water calculations of C ho et alJ (12003) 
suggest that close-in systems like HD 209458b — with strong 
day-night heating contrasts and modest rotation rates — 
occupy a dynamically interesting regime. 

In this article, we report on the results from a multi-layer 
global atmospheric dynamics model of HD 209458b. Our 
results — the existence of a fast superrotating equatorial jet at 
the photosphere that blows the hottest regions downwind — 
qualitatively agree with previous three- dimensional numeri- 
cal simulations by S howman & Guillol ll2002l) . The simula- 
tion presented here, however, adopts more realistic radiative- 
equilibrium temperature profiles and timescales, superior res- 
olution, and a domain that exte nds deeper into the interior . 
In particular, the simulations of S howman & Guillol ( 120021) 
could not accurately predict the day-night temperature differ- 
ence at the photosphere. We here predict this temperature dif- 
ference to b e ~500K, in agreement wit h order of magnitude 
estimates by Showman & Guillot (2002). 

Our calculations have implications for the planet's infrared 
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(IR) light curve. With the first infrared detections of t he tran- 
siting EGPs TrES- 1 and HD 209458b ( Char bonneau. P.. et alJ 
l2005t iDeming et alJl2005T) . observational constraints on the 
meteorologies of close-in giant planets will likely be possible 
over the next two years. 

2. MODEL 

Our model of the general circulation of HP 209458b 
integrates the primitive equations of dynamical meteorol- 
ogy using Version 2 of the ARIES/GEOS Pynamical Core 
( Suarez & Takacs 1995), which we hereafter abbreviate as 
AGPC2. The primitive equations filter vertically propagat- 
ing sound waves but retain horizontally propagating external 
sound waves (called Lamb waves) (Kalnav 2003). 

The grid spacing is 5° x 4° in longitude and latitude, re- 
spectively (^7000 km near the equator). We use 40 vertical 
levels spaced evenly in log-pressure between 1 mbar and 3 
kbar. This spacing implies that we resolve each pressure scale 
height with ~ 2 model layers. The scale height ranges from 
500-1500 k m over the domain of integr ation. 

Foll owing iGuillot et ail ( 119961) and IShowman & Guillot! 
(2002), we assume synchronous rotation. The acceleration 
of gravity, which does not vary significantly over the pres- 
sure range considered, is set to g = 9.42 m s" 2 . We take 
the mean molecular weight and heat capacity to be constant: 
c p = 1.43 x 10 4 Jkg-' K" 1 and n= 1.81 x lO^kgrnol" 1 , which 
neglects the ~ 30% variations in these parameters caused by 
dissociation of molecular hydrogen at the deepest pressures in 
the model. 

The intense stellar irradiation extends the radiative zone of 
HP 209458b all the way down to ~ 1 kbar pressure, with a 
depth ~5-10% of th e planetary radius ( Bur rows et alJl2003l 
Chabrier£La!l2004). Our integrations do not solve the equa- 
tion of radiative transfer directly. Rather, we treat the effects 
of the strong stellar insolation using a Newtonian radiative 
scheme in which the thermodynamic heating rate q [W kg" 1 ] 
is given by 

q _ T(X,cj),p,t)-T eq (X,(j),p) 
c P r md (p) 

which relaxes the model temperature T toward a prescribed 
radiative-equilibrium temperature r eq . In Equation^ A, <j>, p, 



2 



Cooper & Showman 



[ic] 



ST 

■o 

U 
D 



I 



"I - 



, \ j 
- \ i 



I 



■ I 



\ \ "V ~*.-|H 

■ ■ I ■ .... I . 




-ISO -90 90 

Longitude [deg] 



J'.:-. 



6fi = 



ire 



p 1 (Tt 



v a - 

a 

si _ 

3 T 1 



— ^ 



-ISO -90 90 

Longitude [deg] 

— i — ■ — r — | — w r • i i — i — • • ■ — [ — r 




"O 
3 



1 



3 ] 



-ISO 



-90 90 

Longitude [deg] 



I 



2108 
1917 

1725 



FIG. 1 . — Snapshot at 5000 Earth days of simulated temperature (grayscale) 
and winds (arrows) on three isobars: 2.5 mbar, 220 mbar, and 19.6 bars in 
(a), (b), and (c), respectively. The substellar point is at (0,0) in longitude 
and latitude. Peak winds are 9.2, 4.1, and 2.8 km sec -1 in (a), (b), and (c). 
The simulated temperature difference of ~ 500 K at 220 mbar is less than 
the assumed temperature difference of 920 K in radiative-equilibrium due to 
advection of hot material from the dayside to the nightside by eastward winds 
of ~4 km s on this layer. 



and t are longitude, latitude, pressure, and time. The radiative- 
equilibrium temperature T eq and the timescale for relaxation 
to radiative equilibrium T rac j are inputs of the AGDC2. 
We re ly on the radiative-equilibrium calculati ons of 
n (120051) to specify T eq (X,(/),p) and T md (p). Ilro et alJ 
use the multi-wavelength atmosphere code of 
(2000) to calculate the radiative- 



llroet alJ 
(120051) 

Goukenle uque et all 
equilibrium temperature structure of HD 209458b for a single 
vertical column, hereafter denoted as Ti m (p). Ilro et al.l (120051) 
assume globally averaged insolation conditions (i.e., they re- 
distribute the incident solar flux over the entire globe). Their 
calculation includes the opa cities appropriate for a s olar- 
abundance distribution of gas (Anders & Grevesse 1989), in- 
cluding the neutra l alkali metals Na and K. 

Ilro et al 1 (120051) do not consider condensation or the scatter- 
ing and absorption of radiation by silicate clouds, which can 
conce ivably form near the photosphere (see e.g.. lFortnev et al.l 
120031) . Condensates can potentially have a significant effect 
on the planet's radiation balance, depending on the depth at 



which they fo rm, the particle size s, and the vertical extent of 
cloud layers ( Coop eret al.ll2003T) . The net direction of this 
effect (i.e., to warm the atmosphere or to cool it) is as yet 
u nclear and remai ns a subject for future work. 

Ilro et al.l d2005l) compute radiative-relaxation timescales as 
a function of p from 0.01 mbar down to 10 bar by applying 
a Gaussian perturbation to the radiative-equilibrium temper- 
ature profile at each vertical level. We use their radiative- 
relaxation timescales for r md {p) in Equation At pressures 
exceeding 10 bars, radiative relaxation is negligible compared 
to the dynamical timescales considered here. We simply as- 
sume q = on all layers from 10 bar to 3 kbar. 

To account for the longitude-latitude dependence of T eq in 
Equation^ we use a simple prescription. We choose the sub- 
stellar point to be at (A, 0) = (0, 0). On the dayside, we set 

Te q (\ <t>,P) = ^ n1ght,eq(^') + [ ^ st,eq(^')- ^ n'ight■eq(/ :, )] ' COS(A) COS(</>), 

(2) 

where T ss , eq (p) and T n i s h tfiq (p) are the radiative-equilibrium 
temperature profiles of the substellar point and night side (as- 
sumed to be uniform over the dark hemisphere), respectively. 
Equation |2 implies that the hottest T eq profile is at the sub- 
stellar point. On the nightside, we set T eq (X,(f>,p) equal to 

7night,eq(/?)- 

We treat the radiative-equilibrium temperature difference 
between the substellar point and the nightside, AT eq (p) = 
7ss,eq(f) — Tnight.eqO?), as a f ree parameter that is a specified 
functi on of pressure. To determine T sseq (p) and 7riight,eq(f>) 
from Ilro et alJ ( 120051) ' s profile and our specified Ar eq , we 
horizontally average T* q on the top layer of our model over 
the sphere and set it equal to 7£ at that pressure. Based on 
the ~ 1 000 K day-night temperature differences from I ro et alJ 
(2005), we use Ar eq = 1000 K for pressures less than 100 mb 
and decrease it logarithmically with pressure down to 530 K 
at the base of the heated region (10 bar). Newtonian cool- 
ing is a crude approximation to the true radiative transfer, but 
the scheme is computationally fast — hence allowing extensive 
explorations of parameter space — and gives us direct control 
over the model's diabatic heating. 

The model's initial temperature was set to T n i g h teq (p) every- 
where over the globe; there were no initial winds. We set the 
time step equal to 50 s, which is much smaller than the time 
step required for numerical stability according to the Courant- 
Friedrichs-Lewy (CFL) criterion (Kalnav 2003). 

3. RESULTS & DISCUSSION 

By 5000 days of simulation time, the simulation has 
reached a statistical steady-state, at least down to the 3 bar 
level, which is the level above wh ich 99% of the stellar pho- 
tons are absorbed ( Ilro et alJl2005t) . Deeper than 3 bars, the 
kinetic energy continues to increase with time as these layers 
respond to the intense irradiation on relatively long timescales 

Trad ~ 1 >T. 

At pressures less than 10 bars, the model rapidly develops 
strong winds and temperature variability in response to the 
imposed day-night heating contrast. The upper atmosphere 
is nearly in radiative equilibrium (Figure ^i), with tempera- 
ture contrasts of ~ 1000K. This results from the fact that the 
radiative-equilibrium time constant r ra d at 2 mbar is only ~ 1 
hour, which is much shorter than the timescale for winds to 
advect heat across a hemisphere. 

At 2 mbar pressures, supersonic winds exceeding 9 km s -1 
appear at high latitudes, with strong north-south as well as 
east-west flow. Supersonic winds are plausible in the dynami- 
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FIG. 2. — Views of the predicted flux emitted from HD 209458b as it would appear from Earth at four orbital phases: (a) in transit, (b) one-quarter after the 
transit, (c) in secondary eclipse, and (d) one-quarter before the next transit. The planetary rotation axes are vertical, with the superrotating jet seen in FigurefTI 
going from left to right in each panel. The temperature on this layer ranges from 1011 K (dark) to 1526 K (bright). In radiative equilibrium, panel (a) would be 
darkest and (c) would be brightest. The globes show that a difference in observed flux from the planet between the leading and trailing phases of its orbit — panels 
(b) and (d) — is a signature of winds. 

cal regime of HD 209458b due to the immense radiative forc- 
ing from the parent star. For comparison, the planet Neptune 
has high-velocity zona l jets with win d speeds approaching the 
local speed of sound (Limav e~& Sromovskvlll99"il) . but the 
radiative forcing is a million times weaker than it is for HD 
209458b. 

In contrast, the flow near the photosphere near 220 mbar 
is dominated by an eastward jet extending from the equator 
to mid-latitudes (Figure^). The temperature contrasts reach 
~ 500 K at this level, with the hottest part of the atmosphere 
advected ~ 60° downstream from the substellar point by the 
4 km s" 1 eastward jet. Here, the time constant for relaxation 
to the radiative equilibrium r ra d is several Earth days. The 
downstream advection of the hottest regions results from the 
fact that the radiative and advection timescales are compara- 
ble at this level. The small-scale bar-like features visible in 
the middle frame of Figure ^ propagate nearly horizontally 
westward relative to the flow at ^3 km s" 1 , which is close to 
the speed of sound. These features are most consistent with 
Lamb waves (Kalnav 2003, p. 42). 

At 19.6 bar, the equatorial winds remain extremely fast — 
2.8 km s" 1 — but the temperature structure exhibits little lon- 
gitudinal variability (Figure QJ). Our model's circulation re- 
sults entirely from radiative heating occuring at pressures less 
than 10 bars. Therefore, any winds or temperature variabil- 
ity that develops at pressures exceeding 10 bars results solely 
from downward transport of energy from the overlying heated 
layers by vertical advection or wave transport. The layers 
between 3-30 bars contain ~ 70% of the atmosphere's ki- 
netic energy. Lower pressures have fast winds but little mass, 
while the winds drop rapidly to zero at pressures exceeding 
^30 bars. 

The essential assumpti on of the primitive equations is hy- 
drostatic balance l lHoltonl 1 992t lKalnavl2003|) . an approxima- 
tion generally valid for shallow flow. This applies to the case 
of HD 209458b's radiative zone, which has a horizontal to 
vertical aspect ratio of ~ 100. Analysis of the terms in the 
full Navier-Stokes vertical momentum equation reveals that, 
even in the supersonic flow regime obtained here, the verti- 
cal acceleration and curvature terms would have magnitudes 
only ~ 1% and ~ 10% that of the hydrostatic terms at the pho- 
tosphere and top of the model, respectively. The caveat to 
this result is that vertical accelerations can still be important 
for sub-grid scale structures. These simulations also do not 
include the possible effects of vertically propagating shocks, 
which can conceivably dissipate atmospheric kinetic energy. 

To confirm the effects shown in Figure^ we have run addi- 
tional simulations using ID radiative-equilibrium temperature 



profiles from B urrows et"al]j2003l) and lChabrier et alJ <2004l) . 
which are significantly hotter than that o fllro et all d2005l) due 
to differing assumptions about the heat redistribution. We 
have also experimented with changing the value of the free 
parameter AT eq , which controls the strength of the radiative 
forcing in our Newtonian cooling scheme (Equation[Q. We 
have run simulations at Ar eq of 100 K, 250 K, 500 K, 750 K, 
and 1000 K. The essential features of the simulation presented 
here are representative of the results of the other simulations 
performed: they all develop a stable superrotating jet at the 
equator extending to the mid-latitudes and a hot region of at- 
mosphere downwind of the substellar point. But for a given 
Ar eq , the simulated temperature differences are not affected 
by which group's radiative-equilibrium temperature profile is 
used. For values of Ar eq between 500-1000 K, the simu- 
lations all produce temperature contrasts at the photosphere 
ranging from 300-600 K, with peak equatorial winds speeds 
in the range 2-5 km s" 1 . The simulations e mpl oying the atmo- 
spheri c profiles of iBurrows et alJ J2003I) and IChabrier et alJ 
(2004) do, however, produce hotter (by ~500 K) mean pho- 
tospheric temperatures than our nominal case. 

We have also run the model with a different initial con- 
dition. We started this alternate simulation with iden- 
tical input parameters to the simulation shown in Fig- 
ure ^ but with a retrograde zonal wind profile: u = 
-3 km s" 1 cos 4 (0) tan _1 (2 bar//?). We set the meridional wind 
v to zero and set the p ressures to be in gradient-wind balance 
with the initial winds JHoltonll 19921) . After 5000 days of in- 
tegration time, a strong superrotating jet at the equator at 220 
mbar develops in the simulation, with maximum wind speeds 
of3.4kms _1 and temperature contrasts of ^430 K. The simi- 
larity with the simulation presented in Figure^shows that the 
overall flow geometry is not strongly sensitive to the initial 
conditions. 

Our results diffe r from the one-layer shallow-water simula- 
tions of Ch o et alJ (2003), who find that the mean-equatorial 
flow for HD 209458b is westward. Shallow-water turbulence 
simulations consistently produce westward equatorial flow, 
even for planets such as Jupiter and Saturn whose equatorial 
jets are eastward (| Cho & Polvanil9 96llIacono et alJl9 99a b; 
iPeltier & StuhndhoOlt iShowmanl 120041) . This effect mav 
result from the exclusion of three-dimensional momentum- 
transport processes in one-layer models. For example, three- 
dimensional effects are important in allowing equatorial su- 
perrotation in nume rical models of V enus, Earth, and Jupiter 
ilDel Genio & Zhoulfi996l Is aravana n1 H993t lWmiamsl l2003). 
Nevertheless, the shallow-water turbulence models success- 
fully capture essential aspects of mid-latitude jets for giant 
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FIG. 3. — Synthetic lightcurve of the 220 mbar layer of our model atmo- 
sphere for HD 209458b. This simulation predicts peak emission from the 
planet 14 hours before the time of the secondary eclipse. 

planets in our solar system iCho & Polvanill 996). 

4. ORBITAL PHASES AND PREDICTED LIGHTCURVE 

To illustrate the observable effect s of the circu lation, Figure 
|5]shows orthographic projections iSnvdeJl98 7) of the black- 
body flux of HD 209458b on the 220-mbar level. The key 
feature is the hot region downstream from the substellar point, 
which faces Earth after the transit and before the secondary 
eclipse. This pattern differs drastically from that expected in 
the absence of winds, in which case the hottest regions would 
be at the substellar point. 

We show the infrared light curve of the planet as predicted 
by our simulations in Figure[3] The circles represent the total 
flux received at Earth from the planet every 10° in its orbit. 
The points are derived by integration of the radiation inten- 
sity emitted into a solid angle projected in the direction to- 
ward Earth. We assume in this calculation that each column 
of the atmosphere at the photosphere emits radiation with the 
intensity of a blackbody, <tT 4 / it, with temperatures varying as 
shown in Figure^b). We use 223 mbar for the photospheric 
pressure, which is the laye r closest to where the effective tem- 
perature in Iro et al. (2005)'s model equals their computed ac- 
tual temperature. 



The model predicts a phase lead of 60° — or 14 hours — 
between peak radiation from the planet and the secondary 
eclipse, when the illuminated hemisphere faces Earth. Based 
on the ^500 K temperature contrasts shown in Figure |3j we 
derive a ratio of 2.2 between the maximum and minimum flux. 
To confirm the validity of this prediction, we performed anal- 
ogous light curve calculations assuming the IR flux is emit- 
ted from a single pressure level ranging from 150-450 mbar, 
where the photosphere conceivably could be. The magnitude 
of the flux ratio ranges from 3.2 to 1 .4 over this range of pres- 
sures. Our predicted phase lead of 60° shifts by ±20° over 
this range of pressures. These effects are inversely correlated: 
if the photosphere pressure is less than 220 mbar, then the flux 
ratio increases but the phase shift decreases and vice-versa. 

Nevertheless, the possible formation of clouds high in the 
atmosphere of HD 209458b (p < 200 mbar) — which we have 
not treated here — may significantly alter the radiation bud- 
get. Additionally, uncertainties in the input parameters and 
the Newtonian cooling approximation limit the precision with 
which we can predict the magnitude of the effects described 
above. A major future advance in the characterization of 
close-in EGP atmospheres will come from coupling dynam- 
ics and radiative transfer directly. We note, however, that the 
inputs to such a model are likely to be somewhat hypothetical 
as well, given the limited constraints available from observa- 
tions. 
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